The goal of this practical is to introduce you to fitting frailty models in R. We will try out functions from the survival and frailtyEM packages. We will not focus on covariate selection or cover the medical implications of the data analysis. Make sure that you have an up-to-date R installation.
The data that we will use comes from a placebo controlled trial of gamma interferon in chronic granulotomous disease (CGD). The main study based on these data was published as Gallin, J. I., et al. A controlled trial of interferon gamma to prevent infection in chronic granulomatous disease, in The New England Journal of Medicine 324.8 (1991): 509-516. The data set is available in numerous R packages, including the survival package.
The structure of the data is as follows:
id: subject idcenter: the centers where the trial took placerandom: randomization datetreat: treatment arm: rIFN-g or placebosex: female or maleage: age in years at time of study entryheight: height in cm at time of study entryweight: weight in kg at time of study entryinherit: pattern of inheritance: autosomal (autosomal recessive) or X-linkedpropylac: use of prophylactic antibiotics at time of study entry (0 for no, and 1 for yes)hos.cat: institution category: US:other, US:NIH, Europe: Amsterdam and Europe:othertstart, tstop, status: event times from randomization (time 0) to infectionsenum the number of the eventFirst we load the data:
library(survival)
data(cgd)
This loads two data sets: cgd and cgd0. We will use cgd, which is already in the Andersen-Gill format that we discussed in the course. You can find a brief description of this data set in the survival package:
?cgd
head(cgd)
(Q1) Which event is the outcome of interest here? How many individuals are there in the data set? How many individuals are within each center?
length(unique(cgd$id))
## [1] 128
tapply(cgd[cgd$enum == 1,"id"], cgd[cgd$enum == 1, "center"], length)
## Harvard Medical Sch Scripps Institute Copenhagen
## 4 16 4
## NIH L.A. Children's Hosp Mott Children's Hosp
## 26 8 9
## Univ. of Utah Univ. of Washington Univ. of Minnesota
## 4 4 6
## Univ. of Zurich Texas Children's Hosp Amsterdam
## 16 8 19
## Mt. Sinai Medical Ctr
## 4
The outcome of interest is the occurence of recurrent infections after the interventions. It is determined by 3 columns: tstart, tstop and status. We can see that there are 128 individuals in the data set. In terms of centers we see that there are a few large ones (NIH, Amsterdam, Scripps and Zurich) and the rest are smaller.
(Q2) How many events are there in the whole data set?
sum(cgd$status)
## [1] 76
(A2) There are 76 events in the whole data set.
(Q3) How many events are there for each individual? Make a histogram of this. What percentage of individuals experience no serious infections? What are the maximum number of observed events / individual?
hist(tapply(cgd$status, cgd$id, sum), main = "events / individual", xlab = "events")
mean(tapply(cgd$status, cgd$id, sum) == 0)
## [1] 0.65625
(A3) In the histogram (barplot) we see that the most individuals actually have 0 events, that is around 65% of patients. The maximum number of events for an individual is 7.
With recurrent events, often the cumulative intensity is more interesting than the survival. From counting process theory, a counting process \(N(t)\) with intensity \(\lambda(t)\) may be decomposed into: \[ N(t) = \Lambda(t) + M(t). \] This classical result is known as the Doobs-Meyer decomposition and it is fundamental in applying martingale theory to survival analysis. The interpretation here is that \(N(t)\) is the number of events of an individual up to time \(t\), \(M(t)\) is the martingale residual at time \(t\) and \(\Lambda(t)\) is the cumulative intensity (hazard), which for a proportional intensity (hazard) would be \[ \Lambda(t) = \int_0^t \lambda_i(s)ds = \int_0^t \exp(\beta^\top x_i) \lambda_0(s) ds. \] We will not dwelve on this here, but \(M(t)\) has expectation 0, which implies that \(\Lambda(t)\) is equal to the expected number of events of an individual up to time \(t\). When the outcome is a single event (such as death), then \(N(t) \leq 1\) and the intensity \(\lambda(t)\) is the hazard function.
The survival, in the recurrent event case, would be \[ S(t) = \exp(-\Lambda(t)) \] which is not a quantity that can be easily interpreted.
First, we will look at the time to the first event for each individual.
cgd1 <- cgd[cgd$enum==1,]
plot(survfit(Surv(tstop, status) ~ 1, cgd1))
(Q4) What is the interpretation of this curve?
(A4) This is a Kaplan Meier curve for the time to the first event. As we saw before, most individuals actually never get an event. After 200 days, around 80% of individuals still have not had an event.
(Q5) Which variables influence the time to the first event, and how? Tip: you can call summary() on an object for an easier to read output.
coxph(Surv(tstop, status) ~ frailty(id), cgd1)
coxph(Surv(tstop, status) ~ sex + treat + age, cgd1)
coxph(Surv(tstop, status) ~ sex + treat + age + frailty(id), cgd1)
mod_univ <- coxph(Surv(tstop, status) ~ sex + treat + age + inherit + steroids +
propylac + frailty(id), cgd1)
summary(mod_univ)
## Call:
## coxph(formula = Surv(tstop, status) ~ sex + treat + age + inherit +
## steroids + propylac + frailty(id), data = cgd1)
##
## n= 128, number of events= 44
##
## coef se(coef) se2 Chisq DF p
## sexfemale -0.49607 0.50812 0.50812 0.95 1 0.3300
## treatrIFN-g -1.12082 0.34439 0.34439 10.59 1 0.0011
## age -0.03743 0.01851 0.01851 4.09 1 0.0430
## inheritautosomal 0.49277 0.39981 0.39981 1.52 1 0.2200
## steroids 1.10624 0.75469 0.75469 2.15 1 0.1400
## propylac -0.40837 0.42269 0.42269 0.93 1 0.3300
## frailty(id) 0.00 0 0.9400
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexfemale 0.6089 1.6422 0.2249 1.6484
## treatrIFN-g 0.3260 3.0674 0.1660 0.6403
## age 0.9633 1.0381 0.9289 0.9988
## inheritautosomal 1.6368 0.6109 0.7476 3.5837
## steroids 3.0230 0.3308 0.6887 13.2687
## propylac 0.6647 1.5044 0.2903 1.5221
##
## Iterations: 6 outer, 22 Newton-Raphson
## Variance of random effect= 5e-07 I-likelihood = -184.9
## Degrees of freedom for terms= 1 1 1 1 1 1 0
## Concordance= 0.681 (se = 0.047 )
## Likelihood ratio test= 18.34 on 6 df, p=0.005444
(A5) In the last model we can see that treatment and age are the only significant variables. Treatment and age both decrease the time to the first event. Treatment has a very significant effect: it decreases the hazard by about 68%!
(Q6) Do you think that the frailty() statement here is useful? Try to fit the same models without the frailty() statement. Does anything change?
(A6) Probably it will not make any difference: the estimated frailty variance is virtually 0. This means that the fitted model is basically a Cox model. You can see that the results are completely identical if we fit the same model without frailty.
mod_univ_nofr <- coxph(Surv(tstop, status) ~ sex + treat + age + inherit + steroids +
propylac, cgd1)
summary(mod_univ_nofr)
## Call:
## coxph(formula = Surv(tstop, status) ~ sex + treat + age + inherit +
## steroids + propylac, data = cgd1)
##
## n= 128, number of events= 44
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sexfemale -0.49607 0.60892 0.50812 -0.976 0.32893
## treatrIFN-g -1.12082 0.32601 0.34439 -3.255 0.00114 **
## age -0.03743 0.96326 0.01851 -2.022 0.04315 *
## inheritautosomal 0.49277 1.63684 0.39981 1.232 0.21776
## steroids 1.10624 3.02296 0.75469 1.466 0.14270
## propylac -0.40837 0.66473 0.42269 -0.966 0.33398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexfemale 0.6089 1.6422 0.2249 1.6484
## treatrIFN-g 0.3260 3.0674 0.1660 0.6403
## age 0.9633 1.0381 0.9289 0.9988
## inheritautosomal 1.6368 0.6109 0.7476 3.5837
## steroids 3.0230 0.3308 0.6887 13.2687
## propylac 0.6647 1.5044 0.2903 1.5221
##
## Concordance= 0.677 (se = 0.047 )
## Rsquare= 0.133 (max possible= 0.952 )
## Likelihood ratio test= 18.34 on 6 df, p=0.005444
## Wald test = 16.97 on 6 df, p=0.009404
## Score (logrank) test = 18.16 on 6 df, p=0.005839
(Q7) From the output of mod_univ, what do you think that the empirical Bayes estimates of the fraily will look like? Try to find the estimates. You can see the fields corresponding to an object in R by calling str() on that object. On which scale do you think that the frailty estimates are on?
str(mod_univ)
(A7) The frailty estimates are in the $frail field. The estimates are virtually 0, but we can see that a number of them are negative. In fact, they are the logarithm of the empirical Bayes estimates (not the estimates of the log-frailty!).
hist(mod_univ$frail)
(Q8) Do you think that non-proportional covariate effect might be a big problem for this model? Call cox.zph() on it and check the results.
(A8) We can’t really know that unless we test it.
cox.zph(mod_univ)
## rho chisq p
## sexfemale -0.0186 0.0146 0.904
## treatrIFN-g -0.0194 0.0193 0.889
## age 0.1044 0.5180 0.472
## inheritautosomal -0.0168 0.0157 0.900
## steroids 0.0272 0.0273 0.869
## propylac -0.0271 0.0331 0.856
## GLOBAL NA 0.7096 0.994
Nothing appears to be significant which shows that there is no evidence of non-proportionality. If we would have had a significant frailty, probably that there would be some non-proportionality in the Cox model (since the univariate frailty model is confounded with non-proportional covariate effects).
(Q8b) Now back to analyzing the whole data set. In what time scale are the event times measured in this data set? Do you think that the calendar time or the gap time are more relevant in this case?
(A8b) I noticed that I put question 8 twice. Sorry! The time scale is time since intervention here. Calendar time is probably more relevant since there is a clear time origin.
We will start with the marginal models with working independence. Below is the fit using the calendar time scale.
(Q9) What is the left hand side in the argument of coxph? Why was it all right to use only (tstop, status) before?
mod_cal_wi <- coxph(Surv(tstart, tstop, status) ~ sex + treat + age + inherit +
steroids + propylac + cluster(id),
ties = "breslow", cgd)
summary(mod_cal_wi)
## Call:
## coxph(formula = Surv(tstart, tstop, status) ~ sex + treat + age +
## inherit + steroids + propylac + cluster(id), data = cgd,
## ties = "breslow")
##
## n= 203, number of events= 76
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sexfemale -0.73418 0.47990 0.39036 0.43343 -1.694 0.090290 .
## treatrIFN-g -1.05399 0.34855 0.26504 0.30914 -3.409 0.000651 ***
## age -0.04235 0.95853 0.01395 0.01479 -2.863 0.004194 **
## inheritautosomal 0.68068 1.97522 0.27413 0.38858 1.752 0.079820 .
## steroids 1.39243 4.02463 0.56298 0.63176 2.204 0.027521 *
## propylac -0.54460 0.58007 0.30357 0.35074 -1.553 0.120494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexfemale 0.4799 2.0838 0.2052 1.1222
## treatrIFN-g 0.3485 2.8691 0.1902 0.6388
## age 0.9585 1.0433 0.9311 0.9867
## inheritautosomal 1.9752 0.5063 0.9223 4.2303
## steroids 4.0246 0.2485 1.1667 13.8829
## propylac 0.5801 1.7239 0.2917 1.1535
##
## Concordance= 0.702 (se = 0.035 )
## Rsquare= 0.168 (max possible= 0.966 )
## Likelihood ratio test= 37.37 on 6 df, p=1.495e-06
## Wald test = 26.9 on 6 df, p=0.0001515
## Score (logrank) test = 36.78 on 6 df, p=1.94e-06, Robust = 12.62 p=0.0495
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).
(A9) The left hand side is now Surv(tstart, tstop, status). Before, when looking at the first event only, tstart was always 0. That is what R considers if you do not include tstart.
(Q10) We used + cluster(id) to let coxph know which observations to take together. What happens if we would not add that statement? Fit the same model without that statment. Do you notice any differences?
(A10) The difference is that +cluster() gives the robust standard errors (the robust se column in the output), which is a correct estimate when the observations are not independent.
(Q11) Suppose now that we want to go in the other direction - and try a fixed effects model. We can do that by adding instead of + cluster(id), + as.factor(id). Fit this model and examine the output. Do you see anything problematic? How could you explain that?
mod_cal_fe <- coxph(Surv(tstart, tstop, status) ~ sex + treat + age + inherit +
steroids + propylac + as.factor(id),
ties = "breslow", cgd)
mod_cal_fe
## Call:
## coxph(formula = Surv(tstart, tstop, status) ~ sex + treat + age +
## inherit + steroids + propylac + as.factor(id), data = cgd,
## ties = "breslow")
##
## coef exp(coef) se(coef) z p
## sexfemale -6.39e+01 1.69e-28 5.85e+02 -0.11 0.91
## treatrIFN-g -9.59e+01 2.26e-42 7.87e+02 -0.12 0.90
## age 2.91e+00 1.83e+01 2.61e+01 0.11 0.91
## inheritautosomal 2.92e+00 1.85e+01 1.18e+02 0.02 0.98
## steroids -4.98e+01 2.41e-22 4.43e+02 -0.11 0.91
## propylac -7.81e+01 1.18e-34 6.32e+02 -0.12 0.90
## as.factor(id)2 -8.92e+01 1.86e-39 8.16e+02 -0.11 0.91
## as.factor(id)3 -1.15e+01 1.05e-05 1.42e+02 -0.08 0.94
## as.factor(id)4 8.88e+00 7.21e+03 1.24e+02 0.07 0.94
## as.factor(id)5 -9.33e+01 2.94e-41 8.43e+02 -0.11 0.91
## as.factor(id)6 -1.01e+02 1.04e-44 8.36e+02 -0.12 0.90
## as.factor(id)7 -1.08e+02 8.99e-48 9.72e+02 -0.11 0.91
## as.factor(id)8 2.33e+01 1.38e+10 2.31e+02 0.10 0.92
## as.factor(id)9 -1.26e+02 2.57e-55 1.13e+03 -0.11 0.91
## as.factor(id)10 2.92e+01 4.70e+12 2.79e+02 0.10 0.92
## as.factor(id)11 -4.85e+01 8.91e-22 4.60e+02 -0.11 0.92
## as.factor(id)12 2.80e+01 1.49e+12 2.02e+02 0.14 0.89
## as.factor(id)13 8.88e+00 7.21e+03 1.24e+02 0.07 0.94
## as.factor(id)14 -4.57e+01 1.44e-20 4.35e+02 -0.11 0.92
## as.factor(id)15 2.35e+01 1.65e+10 1.94e+02 0.12 0.90
## as.factor(id)16 -1.31e+02 1.78e-57 1.10e+03 -0.12 0.91
## as.factor(id)17 -8.46e+01 1.80e-37 7.66e+02 -0.11 0.91
## as.factor(id)18 -4.73e+01 2.87e-21 4.35e+02 -0.11 0.91
## as.factor(id)19 8.87e+00 7.13e+03 1.69e+02 0.05 0.96
## as.factor(id)20 -1.32e+02 5.69e-58 1.18e+03 -0.11 0.91
## as.factor(id)21 1.47e+01 2.41e+06 1.62e+02 0.09 0.93
## as.factor(id)22 2.33e+01 1.34e+10 2.32e+02 0.10 0.92
## as.factor(id)23 -6.68e+01 9.54e-30 5.25e+02 -0.13 0.90
## as.factor(id)24 8.72e+00 6.12e+03 1.52e+02 0.06 0.95
## as.factor(id)25 -6.45e+01 9.31e-29 5.86e+02 -0.11 0.91
## as.factor(id)26 -5.52e+01 1.07e-24 4.38e+02 -0.13 0.90
## as.factor(id)27 3.49e+01 1.42e+15 3.30e+02 0.11 0.92
## as.factor(id)28 -7.26e+01 3.04e-32 5.89e+02 -0.12 0.90
## as.factor(id)29 3.70e+01 1.19e+16 2.75e+02 0.13 0.89
## as.factor(id)30 -4.72e+01 3.10e-21 4.35e+02 -0.11 0.91
## as.factor(id)31 -5.82e+01 5.36e-26 5.36e+02 -0.11 0.91
## as.factor(id)32 -1.08e+02 2.03e-47 8.96e+02 -0.12 0.90
## as.factor(id)33 5.84e+00 3.44e+02 1.15e+02 0.05 0.96
## as.factor(id)34 -1.62e+02 3.06e-71 1.32e+03 -0.12 0.90
## as.factor(id)35 -1.22e+02 1.05e-53 1.03e+03 -0.12 0.91
## as.factor(id)36 1.50e+01 3.37e+06 1.01e+02 0.15 0.88
## as.factor(id)37 3.78e+01 2.61e+16 3.54e+02 0.11 0.92
## as.factor(id)38 -5.81e+00 3.01e-03 1.17e+02 -0.05 0.96
## as.factor(id)39 -2.35e+02 8.35e-103 1.96e+03 -0.12 0.90
## as.factor(id)40 2.33e+01 1.30e+10 2.32e+02 0.10 0.92
## as.factor(id)41 2.79e+00 1.63e+01 1.01e+02 0.03 0.98
## as.factor(id)42 2.04e+01 7.12e+08 2.09e+02 0.10 0.92
## as.factor(id)43 -8.97e+01 1.13e-39 7.42e+02 -0.12 0.90
## as.factor(id)44 3.99e+01 2.17e+17 3.00e+02 0.13 0.89
## as.factor(id)45 7.85e+01 1.26e+34 7.11e+02 0.11 0.91
## as.factor(id)46 -1.25e+02 3.18e-55 1.13e+03 -0.11 0.91
## as.factor(id)47 -4.06e+01 2.22e-18 2.92e+02 -0.14 0.89
## as.factor(id)48 -1.25e+02 5.42e-55 1.05e+03 -0.12 0.91
## as.factor(id)49 -6.63e+01 1.55e-29 6.12e+02 -0.11 0.91
## as.factor(id)50 -1.20e+02 9.64e-53 1.08e+03 -0.11 0.91
## as.factor(id)51 6.69e+01 1.11e+29 6.08e+02 0.11 0.91
## as.factor(id)52 -1.01e+02 1.39e-44 9.21e+02 -0.11 0.91
## as.factor(id)53 -1.36e+02 1.30e-59 1.16e+03 -0.12 0.91
## as.factor(id)54 -5.81e+00 2.99e-03 1.18e+02 -0.05 0.96
## as.factor(id)55 1.69e+01 2.20e+07 1.17e+02 0.14 0.88
## as.factor(id)56 -2.32e+01 8.36e-11 2.59e+02 -0.09 0.93
## as.factor(id)57 -8.58e+01 5.30e-38 7.29e+02 -0.12 0.91
## as.factor(id)58 -1.48e+02 4.27e-65 1.26e+03 -0.12 0.91
## as.factor(id)59 2.62e+01 2.28e+11 2.57e+02 0.10 0.92
## as.factor(id)60 7.60e+01 9.67e+32 6.86e+02 0.11 0.91
## as.factor(id)61 -1.49e+02 2.82e-65 1.27e+03 -0.12 0.91
## as.factor(id)62 -5.88e+00 2.80e-03 5.22e+01 -0.11 0.91
## as.factor(id)63 -5.88e+01 2.94e-26 5.36e+02 -0.11 0.91
## as.factor(id)64 -2.40e+01 3.84e-11 3.00e+02 -0.08 0.94
## as.factor(id)65 -2.91e+01 2.33e-13 3.05e+02 -0.10 0.92
## as.factor(id)66 -1.48e+02 3.73e-65 1.26e+03 -0.12 0.91
## as.factor(id)67 1.45e+01 2.01e+06 2.04e+02 0.07 0.94
## as.factor(id)68 -4.04e+01 2.96e-18 2.94e+02 -0.14 0.89
## as.factor(id)69 2.62e+01 2.29e+11 2.57e+02 0.10 0.92
## as.factor(id)70 NA NA 0.00e+00 NA NA
## as.factor(id)71 5.52e+01 9.41e+23 5.07e+02 0.11 0.91
## as.factor(id)72 -1.20e+01 6.05e-06 2.32e+02 -0.05 0.96
## as.factor(id)73 -6.97e+01 5.17e-31 5.64e+02 -0.12 0.90
## as.factor(id)74 7.85e+01 1.20e+34 7.12e+02 0.11 0.91
## as.factor(id)75 -7.56e+01 1.52e-33 6.14e+02 -0.12 0.90
## as.factor(id)76 -8.43e+01 2.49e-37 6.92e+02 -0.12 0.90
## as.factor(id)77 -9.30e+01 4.06e-41 7.69e+02 -0.12 0.90
## as.factor(id)78 -1.33e+02 1.21e-58 1.13e+03 -0.12 0.91
## as.factor(id)79 -1.95e+02 2.66e-85 1.68e+03 -0.12 0.91
## as.factor(id)80 3.78e+01 2.52e+16 3.56e+02 0.11 0.92
## as.factor(id)81 2.32e+01 1.23e+10 2.35e+02 0.10 0.92
## as.factor(id)82 -5.52e+01 1.04e-24 4.40e+02 -0.13 0.90
## as.factor(id)83 -7.56e+01 1.53e-33 6.17e+02 -0.12 0.90
## as.factor(id)84 -8.14e+01 4.54e-36 6.67e+02 -0.12 0.90
## as.factor(id)85 -1.25e+02 7.48e-55 9.71e+02 -0.13 0.90
## as.factor(id)86 4.03e+01 3.19e+17 3.00e+02 0.13 0.89
## as.factor(id)87 -9.87e+01 1.37e-43 8.95e+02 -0.11 0.91
## as.factor(id)88 -7.27e+01 2.79e-32 5.90e+02 -0.12 0.90
## as.factor(id)89 2.03e+01 6.74e+08 2.12e+02 0.10 0.92
## as.factor(id)90 -7.60e+01 1.03e-33 6.89e+02 -0.11 0.91
## as.factor(id)91 -7.85e+01 8.33e-35 6.41e+02 -0.12 0.90
## as.factor(id)92 -1.45e+01 4.88e-07 1.68e+02 -0.09 0.93
## as.factor(id)93 -1.54e+02 1.26e-67 1.31e+03 -0.12 0.91
## as.factor(id)94 2.90e+01 4.11e+12 2.83e+02 0.10 0.92
## as.factor(id)95 2.62e+01 2.29e+11 2.57e+02 0.10 0.92
## as.factor(id)96 1.16e+01 1.10e+05 1.50e+02 0.08 0.94
## as.factor(id)97 2.90e+01 4.11e+12 2.83e+02 0.10 0.92
## as.factor(id)98 -5.46e+01 1.94e-24 5.08e+02 -0.11 0.91
## as.factor(id)99 -8.74e+00 1.61e-04 1.36e+02 -0.06 0.95
## as.factor(id)100 2.90e+01 3.88e+12 2.32e+02 0.12 0.90
## as.factor(id)101 4.07e+01 4.64e+17 3.84e+02 0.11 0.92
## as.factor(id)102 -2.91e+01 2.36e-13 2.82e+02 -0.10 0.92
## as.factor(id)103 -1.42e+02 1.42e-62 1.21e+03 -0.12 0.91
## as.factor(id)104 1.75e+01 3.85e+07 1.87e+02 0.09 0.93
## as.factor(id)105 -2.62e+01 4.30e-12 2.60e+02 -0.10 0.92
## as.factor(id)106 -2.62e+01 4.34e-12 1.69e+02 -0.15 0.88
## as.factor(id)107 1.45e+01 2.01e+06 1.72e+02 0.08 0.93
## as.factor(id)108 -1.28e+02 4.12e-56 1.08e+03 -0.12 0.91
## as.factor(id)109 -1.28e+02 2.94e-56 1.08e+03 -0.12 0.91
## as.factor(id)110 3.21e+01 8.53e+13 2.26e+02 0.14 0.89
## as.factor(id)111 -1.28e+02 2.95e-56 1.08e+03 -0.12 0.91
## as.factor(id)112 -2.32e+01 8.04e-11 1.44e+02 -0.16 0.87
## as.factor(id)113 -1.02e+02 6.65e-45 8.46e+02 -0.12 0.90
## as.factor(id)114 -8.14e+01 4.55e-36 6.66e+02 -0.12 0.90
## as.factor(id)115 -2.03e+02 6.44e-89 1.68e+03 -0.12 0.90
## as.factor(id)116 -1.91e+02 7.19e-84 1.58e+03 -0.12 0.90
## as.factor(id)117 -8.72e+01 1.37e-38 7.18e+02 -0.12 0.90
## as.factor(id)118 3.20e+01 7.64e+13 3.12e+02 0.10 0.92
## as.factor(id)119 -5.68e+01 2.10e-25 5.34e+02 -0.11 0.92
## as.factor(id)121 -5.52e+01 1.06e-24 4.26e+02 -0.13 0.90
## as.factor(id)122 -1.02e+02 6.73e-45 8.47e+02 -0.12 0.90
## as.factor(id)123 -9.86e+01 1.52e-43 8.95e+02 -0.11 0.91
## as.factor(id)125 NA NA 0.00e+00 NA NA
## as.factor(id)131 NA NA 0.00e+00 NA NA
## as.factor(id)132 -2.73e+00 6.55e-02 2.61e+01 -0.10 0.92
## as.factor(id)133 NA NA 0.00e+00 NA NA
## as.factor(id)134 NA NA 0.00e+00 NA NA
## as.factor(id)135 NA NA 0.00e+00 NA NA
##
## Likelihood ratio test=178 on 127 df, p=0.00198
## n= 203, number of events= 76
(A11) The estimates of the fixed effect for individuals are either 0 or very large. This is a probelm of collinearity: if we have fixed individual effects, for example, we can’t tell the effect of sex or other variables that are constant within the individual.
(Q12) We can also try to stratify, so that we let each individual have their own baseline hazard. In this case, we would use + strata(id) instead of + cluster(id). What happens? Any idea why?
try(mod_cal_strat <- coxph(Surv(tstart, tstop, status) ~ sex + treat + age + inherit +
steroids + propylac + strata(id),
ties = "breslow", cgd))
(A12) This won’t work. The error is quite funny (This should never happen. Please contact the author.), but the reason is the same as before: since each individual has its own baseline hazard, we can’t decide on covariates effect if they do not vary within the strata.
(Q13) What do you think is the key for the fixed effects and the stratified models to work? In which case could you use those models?
(A13) We would need covariates to vary within the individual. Strata would be a good idea if there would be many events within each cluster/individual. For instance, when there are large centers with many patients. Fixed effects assume proportionality between clusters, and that might work even when there are not so many events per cluster or individual.
Although the coxph and frailty functions are sometimes enough to get an idea of the model fit, we can get some more results with the frailtyEM package. If you don’t have it already installed, you can do that with install.packages("frailtyEM"). The syntax is similar but a bit different from coxph.
(Q18) How do you identify the clusters with the emfrail() function? What is the default frailty distribution? What frailty distributions can be fitted with this function?
library(frailtyEM)
?emfrail
?emfrail_dist
(A18) In frailtyEM we use +cluster() to tell R which observations belong to the same cluster/individual. This plays a similar role with +frailty() from coxph. The default frailty distribution is the gamma, but also positive stable and PVF can be fitted. The PVF also has a parameter m, which by default is -0.5, defining the inverse Gaussian distribution. With positive m, we would obtain compound Poisson distributions.
(Q19) A gamma frailty fit is very similar to what we have seen with coxph. Do you notice any differences in estimates with the coxph fit?
mod_gam_1 <- emfrail(Surv(tstart, tstop, status) ~ sex + treat + age + propylac + inherit +
steroids + cluster(id), cgd)
summary(mod_gam_1)
## Call:
## emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
## age + propylac + inherit + steroids + cluster(id), data = cgd)
##
## Regression coefficients:
## coef exp(coef) se(coef) adjusted se z
## sexfemale -0.804123 0.447480 0.475733 0.477113 -1.690281
## treatrIFN-g -0.991684 0.370951 0.299325 0.299997 -3.313073
## age -0.041893 0.958973 0.016605 0.016605 -2.522952
## propylac -0.601626 0.547920 0.393676 0.393932 -1.528225
## inheritautosomal 0.674073 1.962214 0.334486 0.334486 2.015253
## steroids 1.452841 4.275245 0.754092 0.754893 1.926611
## p
## sexfemale 0.0910
## treatrIFN-g 0.0009
## age 0.0116
## propylac 0.1265
## inheritautosomal 0.0439
## steroids 0.0540
## Estimated distribution: gamma / left truncation: FALSE
##
## Fit summary:
## Commenges-Andersen test for heterogeneity: p-val 0.0449
## (marginal) no-frailty Log-likelihood: -323.606
## (marginal) Log-likelihood: -321.087
## LRT: 1/2 * pchisq(5.04), p-val 0.0124
##
## Frailty summary:
## theta = 2.019 (1.29) / 95% CI: [0.744, 23.631]
## variance = 0.495 / 95% CI: [0.042, 1.344]
## Kendall's tau: 0.199 / 95% CI: [0.021, 0.402]
## Median concordance: 0.195 / 95% CI: [0.02, 0.406]
## E[log Z]: -0.268 / 95% CI: [-0.805, -0.021]
## Var[log Z]: 0.637 / 95% CI: [0.043, 2.572]
## Confidence intervals based on the likelihood function
(A19) The estimates are identical, although what happens behind the scenes is quite different.
(Q20) You noticed probably that here is more output here. The field LRT also shows a likelihood ratio test. From the output shown here, what models does this test contrast? Do you see any other evidence for the presence of frailty from the output? What is \(\theta\) in this case?
(A20) This time the LRT is a test between the model with frailty and the same model without frailty, but with the same covariates. It is indeed significant. The Commenges-Andersen test is a generally less powerful score test for heterogeneity, which is also siginificant in this case. \(\theta\) is the inverse of the frailty variance.
(Q21) Look into the contents of the emfrail fit with str(). Can you find the empirical Bayes frailty estimates? How are they different from the ones from the coxph() fit?
plot(exp(mod_cal_fr$frail), mod_gam_1$frail, xlab = "coxph exp($frail)", ylab = "emfrail $frail")
abline(0,1)
They are identical, just that with
frailtyEM what is in `\(frail\) are the actual estimates, not the logarithm.
(Q22) We can also calculate the empirical Bayes estimates by hand (only for the gamma frailty). The formula is: \[
\widehat{z}_i = \frac{\theta + N_i}{\theta + \Lambda_i}
\] where \(N_i\) is the number of events observed on individual \(i\), \(\Lambda_i\) is the summed up cumulative intensity from cluster \(i\). What do you think \(\theta\) is? How is it related to the estimated frailty variance? Take a look at the summary of mod_gam_1. Here are the ingredients that we need to calculate that:
exp(mod_gam_1$logtheta)
mod_gam_1$nevents_id
mod_gam_1$residuals$group
(Q23) Calculate the \(\widehat{z}_i\) from this. Does it agree with the empirical Bayes estimates from coxph or from frailtyEM?
(A23) Let’s see:
myzi <- (exp(mod_gam_1$logtheta) + mod_gam_1$nevents_id) /
(exp(mod_gam_1$logtheta) + mod_gam_1$residuals$group)
plot(myzi, mod_gam_1$frail)
abline(0,1)
They are the same as the ones from
frailtyEM (and the exp() of the ones from coxph).
(Q24) Calculate the sample mean and variance of the estimated \(\widehat{z}_i\). Does the sample variance agree with the estimated frailty variance? Should they? (A24) Well:
mean(myzi)
## [1] 0.9999277
var(myzi)
## [,1]
## [1,] 0.09718679
The mean is always 1 because of the parametrizations. Note that the sample variance of the estimates is not the estimate of the variance (and there’s no reason why it should be).
We can make several types of plots from an emfrail fit. I will use those based on the ggplot2 package here. More info on those can be found in ?autoplot.emfrail. For conventional plots, ?plot.emfrail describes what can be obtained. The advantage of the ggplot2 plots is that they are objects in R and can be easily edited after they are produced.
Start with a histogram of the frailty estimates:
autoplot(mod_gam_1, type = "hist")
From the histogram it’s quite difficult to tell which individual has which value of the frailty. We can make a so-called catterpillar plot:
autoplot(mod_gam_1, type = "frail")
However, it is still difficult to tell to which individual which frailty value belongs to, since there are so many of them. We can use an interactive visualization for this: install the
plotly package with install.packages("plotly") and then try this:
library(plotly)
ggplotly(autoplot(mod_gam_1, type = "frail"))
(Q25) Which are the individuals with the highest and lowest frailty estimates?
(A25) Individual 60 has the lowest frailty and individual 2 has the highest frailty. If you have plotly installed, then you can just hover the mouse and see the name of each individual.
(Q26) From the output of mod_gam_1 we can already read a number of estimates that correspond to log-hazard ratios. How would a marginal hazard ratio look like? Let’s see this, between a treated and an untreated female, with age = 24, with baseline values for the other covariates. What do you expect that will happen with the marginal hazard ratio? How drastic do you think that the effect of the frailty will be, given the strength of the evidence for the frailty?
(A26) The marginal hazard ratio would be shrunken towards 1. Remember that the frailty was significant (p = 0.012). But I do not think that the shrinking will be extreme - I would maybe expect that if the p-value of the frailty would be very small. However, the amount of shrinking does not only depend on the strength of the frailty effect, but also on the other estimates. The only way to find out for sure is to plot it!
sex <- rep("female", 2)
treat <- c("placebo", "rIFN-g")
age <- c(24, 24)
propylac = c(0, 0)
inherit <- rep("X-linked", 2)
steroids <- rep(0, 2)
newdata <- data.frame(sex, treat, age, propylac, inherit, steroids)
autoplot(mod_gam_1, type = "hr", newdata = newdata)
(Q27) Now suppose that we want to predict the cumulative intensity (hazard) for a certain individual. There comes a question: which one is more relevant, the conditional one or the marginal one? Which one would we observe at the level of the population? Which one would be more relevant for a patient?
(A27) The marginal one is what we would observe at the level of the population. For health economics for example that would be more relevant. The patient however cares about what will happen to her, so the conditional hazard ratio is in that case more relevant.
(Q28) Let’s say we take the first woman, on placebo, from the previous question. What does this return? Can you tell what are the fields in this result? Hint: check with ?predict.emfrail
head(predict(mod_gam_1, newdata = newdata[1,]))
## time lp cumhaz se_logH cumhaz_l cumhaz_r se_logH_adj
## 1 3.5 -1.809551 0.000000000 0.0000000 0.0000000000 0.00000000 0.0000000
## 2 4.0 -1.809551 0.004522549 1.0891757 0.0005348823 0.03823916 1.0892153
## 3 6.0 -1.809551 0.009045098 0.8284345 0.0017833459 0.04587658 0.8284866
## 4 8.0 -1.809551 0.013567648 0.7208586 0.0033029145 0.05573292 0.7209184
## 5 11.0 -1.809551 0.018090197 0.6605330 0.0049566276 0.06602377 0.6605982
## 6 14.0 -1.809551 0.022612746 0.6215334 0.0066879556 0.07645629 0.6216028
## cumhaz_l_a cumhaz_r_a survival survival_l survival_r survival_l_a
## 1 0.0000000000 0.00000000 1.0000000 1.0000000 1.0000000 1.0000000
## 2 0.0005348408 0.03824213 0.9954877 0.9994653 0.9624827 0.9994653
## 3 0.0017831639 0.04588126 0.9909957 0.9982182 0.9551598 0.9982184
## 4 0.0033025272 0.05573945 0.9865240 0.9967025 0.9457917 0.9967029
## 5 0.0049559934 0.06603221 0.9820724 0.9950556 0.9361086 0.9950563
## 6 0.0066870462 0.07646669 0.9776410 0.9933344 0.9263934 0.9933353
## survival_r_a survival_m survival_m_l survival_m_r survival_m_l_a
## 1 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 2 0.9624799 0.9954927 0.9994653 0.9628270 0.9994654
## 3 0.9551554 0.9910157 0.9982190 0.9556504 0.9982192
## 4 0.9457855 0.9865688 0.9967052 0.9465065 0.9967056
## 5 0.9361007 0.9821516 0.9950617 0.9370983 0.9950623
## 6 0.9263838 0.9777639 0.9933453 0.9277026 0.9933462
## survival_m_r_a cumhaz_m cumhaz_m_l cumhaz_m_r cumhaz_m_l_a
## 1 1.0000000 0.000000000 0.0000000000 0.00000000 0.0000000000
## 2 0.9628242 0.004517491 0.0005348114 0.03788152 0.0005347699
## 3 0.9556461 0.009024896 0.0017825587 0.04536308 0.0017823769
## 4 0.9465004 0.013522259 0.0033002155 0.05497749 0.0032998289
## 5 0.9370906 0.018009626 0.0049505528 0.06496712 0.0049499201
## 6 0.9276933 0.022487041 0.0066769021 0.07504407 0.0066759957
## cumhaz_m_r_a
## 1 0.00000000
## 2 0.03788444
## 3 0.04536766
## 4 0.05498385
## 5 0.06497530
## 6 0.07505409
This returns a prediction of the cumulative hazard and survival (rather, \(\exp(-\Lambda(t))\)), both marginal and conditional, with 95% confidence intervals.
(Q29) We can look at the plot for the cumulative intensity of our patient of interest. How do you expect that the marginal cumulative intensity will compare to the conditional cumulative intensity?
autoplot.emfrail(mod_gam_1, type = "pred", newdata = newdata[1,])
The dotted lines represent a 95% confidence band.
(A29) Usually the marginal hazard/intensity is “dragged down” by the frailty.
(Q30) Let’s play a bit with the distribution argument in emfrail. I’ll fit 3 more distributions here. Which one is the inverse Gaussian? (hint: ?emfrail_dist())
mod_pvf1 <- emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
age + propylac + inherit + steroids + cluster(id),
distribution = emfrail_dist(dist = "pvf"),
data = cgd)
mod_pvf2 <- emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
age + propylac + inherit + steroids + cluster(id),
distribution = emfrail_dist(dist = "pvf", pvfm = 0.5),
data = cgd)
mod_stab <- emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
age + propylac + inherit + steroids + cluster(id),
distribution = emfrail_dist(dist = "stable"),
data = cgd)
(A30) The default in emfrail_dist(dist = "pvf") is to have pvfm = -0.5. That is the inverse Gaussian distribution. The second model has a compound Poisson distribution.
(Q31) Look at the output of the positive stable frailty model. Is the frailty significant here? Compare the estimates with those from the gamma frailty model and the marginal model with + cluster(id). To which ones is it closer? Is there anything extra in the output? What do you think it means?
summary(mod_stab)
## Call:
## emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
## age + propylac + inherit + steroids + cluster(id), data = cgd,
## distribution = emfrail_dist(dist = "stable"))
##
## Regression coefficients:
## coef exp(coef) se(coef) adjusted se z
## sexfemale -0.669860 0.511780 0.436831 0.441698 -1.533453
## treatrIFN-g -1.045188 0.351626 0.293880 0.294395 -3.556511
## age -0.043091 0.957824 0.015464 0.015505 -2.786602
## propylac -0.599076 0.549319 0.333476 0.340483 -1.796458
## inheritautosomal 0.634477 1.886036 0.330350 0.333529 1.920620
## steroids 1.400698 4.058033 0.609885 0.610805 2.296659
## p
## sexfemale 0.1252
## treatrIFN-g 0.0004
## age 0.0053
## propylac 0.0724
## inheritautosomal 0.0548
## steroids 0.0216
## Estimated distribution: stable / left truncation: FALSE
##
## Fit summary:
## Commenges-Andersen test for heterogeneity: p-val 0.0449
## (marginal) no-frailty Log-likelihood: -323.606
## (marginal) Log-likelihood: -323.416
## LRT: 1/2 * pchisq(0.38), p-val 0.269
##
## Frailty summary:
## theta = 24.099 (39.25) / 95% CI: [4.736, Inf]
## Kendall's tau: 0.04 / 95% CI: [0, 0.174]
## Median concordance: 0.038 / 95% CI: [0, 0.171]
## E[log Z]: 0.024 / 95% CI: [0, 0.122]
## Var[log Z]: 0.139 / 95% CI: [0, 0.768]
## Attenuation factor: 0.96 / 95% CI: [0.826, 1]
## Confidence intervals based on the likelihood function
(A31) This shows that in the positive stable (PS) model the frailty is not significant. In general it is a bad idea to fit the PS model when individuals do not have at least 2 events, since this distribution can not be estimated in that case.
(Q32) Optional: how do you think that the marginal and conditional hazard ratios would look like, if we did the same thing with mod_stab? Try it!
(A32) The marginal and conditional hazard ratios would be parallel lines, but quite close to each other, since the frailty was not significant.
(Q33) How about model mod_pvf2? Notice anything in its summary that wasn’t there before?
summary(mod_pvf2)
## Call:
## emfrail(formula = Surv(tstart, tstop, status) ~ sex + treat +
## age + propylac + inherit + steroids + cluster(id), data = cgd,
## distribution = emfrail_dist(dist = "pvf", pvfm = 0.5))
##
## Regression coefficients:
## coef exp(coef) se(coef) adjusted se z
## sexfemale -0.816874 0.441811 0.476519 0.478196 -1.714251
## treatrIFN-g -0.985376 0.373299 0.298845 0.300314 -3.297282
## age -0.042005 0.958865 0.016567 0.016568 -2.535387
## propylac -0.598970 0.549377 0.393811 0.393876 -1.520958
## inheritautosomal 0.681071 1.975994 0.332952 0.332952 2.045553
## steroids 1.467207 4.337103 0.764430 0.765591 1.919347
## p
## sexfemale 0.0865
## treatrIFN-g 0.0010
## age 0.0112
## propylac 0.1283
## inheritautosomal 0.0408
## steroids 0.0549
## Estimated distribution: pvf / left truncation: FALSE
## PVF m = 0.5
##
## Fit summary:
## Commenges-Andersen test for heterogeneity: p-val 0.0449
## (marginal) no-frailty Log-likelihood: -323.606
## (marginal) Log-likelihood: -321.048
## LRT: 1/2 * pchisq(5.12), p-val 0.0118
##
## Frailty summary:
## theta = 2.058 (1.25) / 95% CI: [0.811, 22.007]
## variance = 0.486 / 95% CI: [0.045, 1.233]
## Estimated mass at 0: 0.002081216
## Confidence intervals based on the likelihood function
(A33) The estimated mass at 0 field is new. That is a feature of the compound Poisson distributions.
(Q34) Compare the likelihoods of all the frailty models we fitted here. Which one is the highest? Was this to be expected? Think of one of the first questions, about how many individuals have how many events.
lapply(list(mod_gam_1, mod_stab, mod_pvf1, mod_pvf2), logLik)
## [[1]]
## 'log Lik.' -321.0875 (df=7)
##
## [[2]]
## 'log Lik.' -323.4158 (df=7)
##
## [[3]]
## 'log Lik.' -321.292 (df=7)
##
## [[4]]
## 'log Lik.' -321.0475 (df=7)
(A34) The compound Poisson has the highest log-likelihood out of all the models. That is to be expected in some sense, since we have seen that most individuals have 0 events: so an estimated mass at 0 is to be expected.
Lastly, we will note that the estimated mass at 0 depends heavily of the pvfm parameter, that defines the large classes of distribution from within the PVF family. This is not estimated by frailtyEM, but rather set by the user. We can however try to estimate it, just by trying out a number of values.
(Q35) I collect here the log-likelihood values for every distribution (for different values of pvfm). What is the value of m for which the highest log-likelihood of the model is obtained?
mvals <- seq(from = 0.1, to = 2, by = 0.2)
models <- lapply(mvals, function(x) emfrail(formula = Surv(tstart, tstop, status) ~ sex +
treat + age + propylac + inherit + steroids + cluster(id), data = cgd,
distribution = emfrail_dist(dist = "pvf", pvfm = x)) )
likelihoods <- sapply(models, function(x) x$loglik[2])
plot(mvals, likelihoods)
(A35) An
m around 1 gives the best fit.
Another option with recurrent events data is to analyze it in gap-time. For this we will have to add another column in our data:
cgd$gap <- cgd$tstop - cgd$tstart
The right hand side of the formulas we used before would stay the same, but the left hand side should now be:
mod_gap_fr <- coxph(Surv(gap, status) ~ sex + treat + age + inherit +
steroids + propylac + frailty(id),
ties = "breslow", cgd)
summary(mod_gap_fr)
## Call:
## coxph(formula = Surv(gap, status) ~ sex + treat + age + inherit +
## steroids + propylac + frailty(id), data = cgd, ties = "breslow")
##
## n= 203, number of events= 76
##
## coef se(coef) se2 Chisq DF p
## sexfemale -0.88169 0.53542 0.42001 2.71 1.00 0.1000
## treatrIFN-g -1.05828 0.33091 0.27776 10.23 1.00 0.0014
## age -0.04506 0.01857 0.01487 5.89 1.00 0.0150
## inheritautosomal 0.72385 0.37952 0.27900 3.64 1.00 0.0560
## steroids 1.55948 0.91448 0.62636 2.91 1.00 0.0880
## propylac -0.72039 0.46262 0.33332 2.42 1.00 0.1200
## frailty(id) 57.94 39.03 0.0260
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexfemale 0.4141 2.4150 0.1450 1.1826
## treatrIFN-g 0.3471 2.8814 0.1814 0.6638
## age 0.9559 1.0461 0.9218 0.9914
## inheritautosomal 2.0624 0.4849 0.9802 4.3393
## steroids 4.7563 0.2102 0.7923 28.5546
## propylac 0.4866 2.0552 0.1965 1.2048
##
## Iterations: 6 outer, 31 Newton-Raphson
## Variance of random effect= 0.9728984 I-likelihood = -340.5
## Degrees of freedom for terms= 0.6 0.7 0.6 0.5 0.5 0.5 39.0
## Concordance= 0.852 (se = 0.036 )
## Likelihood ratio test= 119.6 on 42.52 df, p=2.925e-09
(Q36) Why can’t we use Surv(tstop, status)? What would that correspond to? What is the time scale now?
(A36) Using Surv(tstop, status) would be equivalent with saying that each observation is from 0 to the tstop, and that would count too much at-risk time. We are using the time scale time since previous event now.
Look at the first individual:
cgd[cgd$id == 1, c("treat", "sex", "age", "inherit", "gap", "status")]
## treat sex age inherit gap status
## 1 rIFN-g female 12 autosomal 219 1
## 2 rIFN-g female 12 autosomal 154 1
## 3 rIFN-g female 12 autosomal 41 0
In mod_gap_fr, R doesn’t know that we are talking about recurrent events, since the syntax is exactly the same as for clustered data. (Q37) Do you think that matters here?
(A37) It does not actually! Clustered failures and recurrent events in gap-time are actually estimated in the same way.
(Q38) Below is a Kaplan-Meier curve. What do you think that this means here? Do you think that, if we analyze the models in gap time, the survival is more meaningful?
plot(survfit(coxph(Surv(gap, status) ~ 1, cgd)))
(A38) This curve is the “survival” of the gaptime - it tells us for example that about 70% of the gaptimes are large than 200 (we look at time 200 on the x axis.)
(Q39) Compare the estimated conditional hazard ratios (coefficients) between mod_cal_fr and mod_gap_fr. Are there similarities? Are the differences notable from a clinical point of view, do you think? Would you expect to have large differences between the two models? Do you think anything is lost by looking at the gap times instead of the calendar time?
(A39) The clinical conclusions are pretty similar, although the estimates are different and have a different meaning now. For example, if the treatment increases the intensity of recurrent event process in calendar time, we would also expect that it makes the times between events longer. However, there is no 1-1 correspondence in general between these estimates.
(Q40) Furthermore, compare the frailty estimates from the gap time model and the calendar time model. Do you expect them to be similar? What is different between them? Take a look again at the formula for the gamma empirical Bayes estimates that we looked at before.
plot(exp(mod_gap_fr$frail), exp(mod_cal_fr$frail), xlab = "gaptime frailty", ylab = "calendartime frailty")
abline(0,1)
(A40) They are definetely positively correlated. If a high frailty individual has a high rate of events in calendar time, that would also correspond to shorter gaptimes. The numerator of the formula still contains the number of events for the individual, which stays unchanged in the two models.